#include "stdafx.h"
#include <cuda_runtime.h> 
#include <cutil_inline.h> 

//*****************************************************************************
__constant__ double		d_htcInter[NUM_MATERIAL * NUM_MATERIAL];
__constant__ double		d_density[NUM_MATERIAL];
__constant__ double		d_sh[NUM_MATERIAL];
__constant__ double		d_initTemp[NUM_MATERIAL];

//*****************************************************************************
int setConstant(material substance[], double h_htcInter[])
{
	cudaMemcpyToSymbol(d_htcInter, h_htcInter, 
		sizeof(double) * NUM_MATERIAL * NUM_MATERIAL, 0, cudaMemcpyHostToDevice);

	double h_density[NUM_MATERIAL], h_sh[NUM_MATERIAL], h_initTemp[NUM_MATERIAL];
	for(int i = 0; i < NUM_MATERIAL; i++)
	{
		h_density[i] = substance[i].density;
		h_sh[i] = substance[i].sh;
		h_initTemp[i] = substance[i].initTemp;
	}
	cudaMemcpyToSymbol(d_density, h_density, 
		sizeof(double) * NUM_MATERIAL, 0, cudaMemcpyHostToDevice);
	cudaMemcpyToSymbol(d_sh, h_sh, 
		sizeof(double) * NUM_MATERIAL, 0, cudaMemcpyHostToDevice);
	cudaMemcpyToSymbol(d_initTemp, h_initTemp, 
		sizeof(double) * NUM_MATERIAL, 0, cudaMemcpyHostToDevice);
	return 0;
}

//*****************************************************************************
int getGrid(int size)
{
	return (size + BLOCK_SIZE - 1) / BLOCK_SIZE;
}

//*****************************************************************************
__global__ void kernelInitTemperatureField(int size, int *d_grid, 
										   double *d_tempCurr, double *d_tempNext)
{
	int id = blockIdx.x * blockDim.x + threadIdx.x;
	int initTemp;
	if (id < size)
	{
		initTemp = d_initTemp[d_grid[id]];
		d_tempCurr[id] = initTemp;
		d_tempNext[id] = initTemp;
	}
}

int initTemperatureField(int *d_grid, int row, int col, 
						 double *d_tempCurr, double *d_tempNext) 
{
	dim3 block(BLOCK_SIZE);
	dim3 grid(getGrid(row * col));
	kernelInitTemperatureField<<<grid, block>>>
		(row * col, d_grid, d_tempCurr, d_tempNext);

	//4156tempCurr[0][0] = tempNext[0][0] = substance[0].initTemp;

	return 1;
}

//*****************************************************************************
int outputDeviceTemperatureField(double *d_temp, int row, int col, double time)
{
	double		*h_testDouble = (double *)malloc(sizeof(double) * row * col);
	cudaMemcpy(h_testDouble, d_temp, sizeof(double) * row * col, cudaMemcpyDeviceToHost);
	outputTemperatureField(h_testDouble, row, col, time);
	return 0;
}

//*****************************************************************************
__global__ void kernelTimeStep(int col, int col_2, int matrixSize, double dxdx, 
							   int *d_grid, double *d_timeStep)
{
	int id = blockIdx.x * blockDim.x + threadIdx.x;
	if (id < matrixSize)
	{
		int r = id / col_2;
		int c = id - r * col_2;
		int location = (1 + r) * col + (1 + c);
		int k = d_grid[location];
		d_timeStep[id] = d_density[k] * d_sh[k] * dxdx
						/ (d_htcInter[k * NUM_MATERIAL + d_grid[location - 1]] 
						+ d_htcInter[k * NUM_MATERIAL + d_grid[location + 1]]
						+ d_htcInter[k * NUM_MATERIAL + d_grid[location - col]] 
						+ d_htcInter[k * NUM_MATERIAL + d_grid[location + col]]);
	}
}

//kernel of min of an Array
__global__ void kernelGetMin(const int count, double *d_array)
{
	int id = blockIdx.x * blockDim.x + threadIdx.x;
	if (id < count)
		if (d_array[id] > d_array[id + count]) d_array[id] = d_array[id + count];
}

double getArrayMin(double *d_array, int arrayLength)
{
	dim3		block(BLOCK_SIZE);
	dim3		grid(getGrid(arrayLength / 2));
	int			count = (int)(arrayLength / 2);
	while(count > 1)
	{
		kernelGetMin<<<grid, block>>>(count, d_array);
		count = (count + 1) / 2;
	}
	double		first, last;
	cudaMemcpy(&first, d_array, sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(&last, d_array + arrayLength - 1, sizeof(double), cudaMemcpyDeviceToHost);
	
	return first < last ? first : last;
}

double calcTimeStep(int *d_grid, int row, int col, double dx) 
{
	double		*d_timeStep;
	int			matrixSize = (row - 2) * (col - 2);
	cudaMalloc((void**)&d_timeStep, matrixSize * sizeof(double));
	dim3		block(BLOCK_SIZE);
	dim3		grid(getGrid(matrixSize));

	kernelTimeStep<<<grid, block>>>(col, col - 2, matrixSize, dx * dx, d_grid, d_timeStep);

	return getArrayMin(d_timeStep, matrixSize);
}

double fixTimeStep(double maxTimeStep, double factor) 
{
	return maxTimeStep * factor;
}

//*****************************************************************************
__global__ void kernelCalcNextTemp(double timeStep, int col, int col_2, int matrixSize, double dxdx,
								   int *d_grid, double *d_tempCurr, double *d_tempNext)
{
	int id = blockIdx.x * blockDim.x + threadIdx.x;	
	if (id < matrixSize)
	{
		int r = id / col_2;
		int c = id - r * col_2;
		int location = (1 + r) * col + (1 + c);
		int k = d_grid[location];
		double temp = d_tempCurr[location];
		d_tempNext[location] = temp + timeStep / d_density[k] / d_sh[k] / dxdx * 
			(d_htcInter[k * NUM_MATERIAL + d_grid[location - 1]] * (d_tempCurr[location - 1] - temp) +
			d_htcInter[k * NUM_MATERIAL + d_grid[location + 1]] * (d_tempCurr[location + 1] - temp) +
			d_htcInter[k * NUM_MATERIAL + d_grid[location - col]] * (d_tempCurr[location - col] - temp) +
			d_htcInter[k * NUM_MATERIAL + d_grid[location + col]] * (d_tempCurr[location + col] - temp));
	}
}

double fourierIteration(int *d_grid, int row, int col, 
						double timeStep, double dx, terminationCondition2D tc, 
						double *d_tempCurr, double *d_tempNext) 
{
	double		*swap;
	double		dxdx = dx * dx;
	double		cycle = 0;
	int			matrixSize = (row - 2) * (col - 2);
	dim3		block(BLOCK_SIZE);
	dim3		grid(getGrid(matrixSize));
	double		keyPointTemp;
	cudaMemcpy(&keyPointTemp, d_tempCurr + tc.y * col + tc.x, sizeof(double), cudaMemcpyDeviceToHost);

	outputDeviceTemperatureField(d_tempCurr, row, col, cycle);
	printf("center of plate - %f @ time - %f\n", keyPointTemp, cycle);
	while(keyPointTemp > tc.terminationTemp)
	{
		cycle += timeStep;

		kernelCalcNextTemp<<<grid, block>>>(timeStep, col, col - 2, matrixSize, dx * dx, d_grid, d_tempCurr, d_tempNext);
	
		swap = d_tempCurr; d_tempCurr = d_tempNext; d_tempNext = swap;

		cudaMemcpy(&keyPointTemp, d_tempCurr + tc.y * col + tc.x, sizeof(double), cudaMemcpyDeviceToHost);

		if ((int)(cycle / OUTPUT_FREQUENCY) > (int)((cycle - timeStep) / OUTPUT_FREQUENCY)) 
		{
			outputDeviceTemperatureField(d_tempCurr, row, col, cycle);
			printf("center of plate - %f @ time - %f\n", keyPointTemp, cycle);
		}
	}
	//outputTemperatureField(tempCurr, row, col, cycle + OUTPUT_FREQUENCY * 2);
	printf("center of plate - %f @ time - %f\n", keyPointTemp, cycle);
	return cycle;
	
}

//*****************************************************************************
//fourier method realization
double fourierCuda2D(int argc, char* argv[], 
					 material *substance, int *grid, int row, int col,
					 double htcInter[], 
					 double tempEdge,  double dx, terminationCondition2D tc)
{
	 //init cuda device
	if (cutCheckCmdLineFlag(argc, (const char**)argv, "device"))
		cutilDeviceInit(argc, argv);
	else
		cudaSetDevice(cutGetMaxGflopsDeviceId());

	//timer
	LARGE_INTEGER	nStartCounter, nStopCounter;
	QueryPerformanceCounter(&nStartCounter);//timer starts

	setConstant(substance, htcInter);

	int			*d_grid;
	int			matrixSize = row * col;
	cudaMalloc((void**)&d_grid, matrixSize * sizeof(int));
	cudaMemcpy(d_grid, grid, sizeof(int) * row * col, cudaMemcpyHostToDevice);

	double		*d_tempCurr;
	double		*d_tempNext;
	cudaMalloc((void**)&d_tempCurr, matrixSize * sizeof(double));
	cudaMalloc((void**)&d_tempNext, matrixSize * sizeof(double));

	initTemperatureField(d_grid, row, col, d_tempCurr, d_tempNext);
	//outputTemperatureField(tempCurr, row, col, 0);

	QueryPerformanceCounter(&nStopCounter);//timer ends
	double nTime1 = getCostTime(nStartCounter, nStopCounter);

	double maxTimeStep = calcTimeStep(d_grid, row, col, dx);
	printf("max time step = %f \n", maxTimeStep);
	double timeStep = fixTimeStep(maxTimeStep, TIME_STEP_FACTOR);
	printf("time step = %f \n", timeStep);

	QueryPerformanceCounter(&nStopCounter);//timer ends
	double nTime2 = getCostTime(nStartCounter, nStopCounter);

	double cycle = fourierIteration(d_grid, row, col, 
									timeStep, dx, tc,									
									d_tempCurr, d_tempNext);

	QueryPerformanceCounter(&nStopCounter);//timer ends
	double nTime3 = getCostTime(nStartCounter, nStopCounter);

	return cycle;
}